#include <iostream>
#include <time.h>
#include <stdio.h>
#include "Substance.h"
#include "float.h"
#include <stddef.h>

#include "common/common_types.h"

/////////////////////////////////////
const int IN_PRD = 30;

const float maxR = 0.35f;

/////////////////////////////////////

struct Measur {
    float tempTubeAv;
    float tempTubeMA[MAP];

    float tempOutAv;
    float tempOutMA[MAP];

    float tempInAv;
    float tempInMA[MAP];

    float prsInAv;
    float prsInMA[MAP];

    float streamAv;
    float streamC;

    float prsOutAv;
    float prsOutMA[MAP];
};
////////////////////////////////////////////

Molecula** sys;

Cell ***cells;

////////////////////////////////////////////////
inline float getMA(float * maBuf) {
    float res = 0;
    for (int i = 0; i < MAP; i++) {
        res += maBuf[i];
    }
    return res / MAP;
}

void mpiCalc(int rnk, int * shift, mpi_toSend * packet, float tau) __attribute__((noinline));
void mpiCalc(int rnk, int * shift, mpi_toSend * packet, float tau) {
    int s = 0;
    short sh_ing = 0;
    for (int r = 0; r < rnk; r++) {
        while (shift[s] == 0) {
            s++;
            sh_ing = 1;
        }

        shift[s]--;
        int start = r + 1 + sh_ing * shift[s];

        Vector3 p = packet[r].Pos;
        Vector3 v = packet[r].Vel;
        for (int i = start; i < rnk; i++) {
            if (packet[r].sysID == packet[i].sysID) {
                break;
            }

            Vector3 p_c = p - packet[i].Pos;
            const IFr& a = *Particle::ForceArr[packet[r].compId][packet[i].compId];
            float dist = p_c.length();

            if (dist < 2 * a.Potential(0, 'R')) {
                float M1 = packet[r].M;
                float M2 = packet[i].M;
                Vector3 cross_n = p_c / dist;
                float Meff = M1 * M2 / (M1 + M2);

                float accel = Simple_Calc(a.getForceInfo(), Meff, dist, tau);

                Vector3 v_accel = cross_n * accel;

                Vector3 v_accel1 = v_accel * (Meff / M1);
                Vector3 v_accel2 = v_accel * (-Meff / M2);

                packet[r].Accel += v_accel1;
                packet[i].Accel += v_accel2;
            }
        }
    }
}

void Boxing(int Num_Tube, int totalParticles, int X_B, int Y_B, int Z_B, float hX, float hY, float hZ, const Vector3& L_Bound) __attribute__((noinline));
void Boxing(int Num_Tube, int totalParticles, int X_B, int Y_B, int Z_B, float hX, float hY, float hZ, const Vector3& L_Bound) {
    for (int k = 0; k < Z_B + 2; k++) {
        for (int j = 0; j < Y_B + 2; j++) {
            for (int i = 0; i < X_B + 2; i++) {
                cells[k][j][i].curIn = cells[k][j][i].stIn;
            }
        }
    }

    for (int q = Num_Tube; q < totalParticles; q++) {
        if (sys[q] != 0) {
            int molR = sys[q]->getRank();
            for (int z = 0; z < molR; z++) {
                Component* b = sys[q]->getComponent(z);
                Vector3 v = b->getPos() - L_Bound;

                int I = v.x / hX + 1;
                int J = v.y / hY + 1;
                int K = v.z / hZ + 1;

                if (I < (X_B + 2) && J < (Y_B + 2) && K < (Z_B + 2) && I > 0 && J > 0 && K > 0) {
                    if (cells[K][J][I].curIn < (int) MAX_PARTICLES) {
                        cells[K][J][I].inCell[cells[K][J][I].curIn] = b;
                        cells[K][J][I].curIn += 1;
                    }
                } else {
                    sys[q] = 0;
                    break;
                }
            }
        }

    }
}
//////////////////////////////////////////////
int main() {
    srand((float) time(NULL) * (float) time(NULL));

    FILE* Log = fopen("./IO/Log.txt", "w");
    if (!Log) {
        return -2;
    }
    FILE* InfoFile = fopen("./IO/inInfo1.txt", "r");
    if (!InfoFile) {
        return -2;
    }
    FILE* Moleculs = fopen("./IO/Moleculs.txt", "r");
    if (!Moleculs) {
        return -2;
    }
    FILE* InterMol = fopen("./IO/InterMol.txt", "r");
    if (!InterMol) {
        return -2;
    }
    MainInit mIn = MainInit(InfoFile, Moleculs, InterMol, Log);
    mIn.initAbout();
    mIn.initSub();
    FileReader myin(InfoFile);

    Component::InitComponent(mIn.getSub());

    float l_x = mIn.Tube->Len + 2.0f + 2 * maxR;
    float l_y = mIn.Tube->Rad * 2 + 2 * maxR;
    float l_z = l_y;

    int X_B = ceil(l_x / maxR / 2);
    if (X_B > 20) {
        X_B = X_B / 2;
    }

    int Y_B = ceil(l_y / maxR / 2);
    int Z_B = ceil(l_z / maxR / 2);

    cells = new Cell**[Z_B + 2];
    for (int k = 0; k < Z_B + 2; k++) {
        cells[k] = new Cell*[Y_B + 2];
        for (int j = 0; j < Y_B + 2; j++) {
            cells[k][j] = new Cell[X_B + 2];
        }
    }

    Vector3 L_Bound(-2 * maxR, -l_y / 2 - maxR * 2, -l_z / 2 - maxR * 2);
    Vector3 R_Bound(l_x, l_y / 2 + maxR * 2, l_z / 2 + maxR * 2);

    float hX = (R_Bound.x - L_Bound.x) / (X_B - 1);
    float hY = (R_Bound.y - L_Bound.y) / (Y_B - 1);
    float hZ = (R_Bound.z - L_Bound.z) / (Z_B - 1);

    //////////////////////////////////////////////////////////////////////
    const SubListItem& sub = mIn.getSub();
    MolDscr * mDscr = mIn.getMDscr();

    bool ok;

    float prs = myin.readFloat(PRS, &ok);
    if (!ok) {
        fprintf(Log, "cannot read %s\n", PRS);
        return 2;
    }

    float tmk = myin.readFloat(TEMK, &ok);
    if (!ok) {
        fprintf(Log, "cannot read %s\n", TEMK);
        return 2;
    }

    float ro = myin.readFloat(RO, &ok);
    if (!ok) {
        fprintf(Log, "cannot read %s\n", RO);
        return 2;
    }

    Flow1 inPool = Flow1(sub, mDscr, prs, tmk, ro, mIn.Tube->Rad);
    float Sht = inPool.L_TS;
    //////////////////////////////////////////////////////////////////////
    float t_mass = ro * 1e-18 * mIn.Tube->Vol;
    sys = new Molecula*[(int) ceil(t_mass / inPool.getAvMass() * 20)];

    sys[0] = new TMolecula(0, mIn.Tube, mIn.getSub());
    sys[0]->setPos(Vector3(Sht, 0.0, 0), Quaternion(0, 1, 0, 0));
    int moleculesCount = 1;
    int Num_Tube = 1;
    //////////////////////////////////////////////////////////////////////
    pBound * bound = new pBound[4];
    bound[0] = new RBound(Vector3(-0.5, 0, 0), Vector3(1, 0, 0), mIn.Tube->Rad * 1.6f, 10.0f);
    bound[1] = new PlBound(Vector3(-0.5, 0, 0), Vector3(1, 0, 0), 0.0f, mIn.Tube->Rad * 2.0f, 1000);
    bound[2] = new PlBound(Vector3(Sht + 0.2, 0, 0), Vector3(1, 0, 0), mIn.Tube->Rad + 0.15f, mIn.Tube->Rad * 2.0f, 1000);
    bound[3] = new PlBound(Vector3(Sht - 0.2 + mIn.Tube->Len, 0, 0), Vector3(1, 0, 0), mIn.Tube->Rad + 0.15f, mIn.Tube->Rad * 2.0f, 1000);
    ///////////////////////////////////////////////////////////////////////
    for (int k = 0; k < Z_B + 2; k++) {
        for (int j = 0; j < Y_B + 2; j++) {
            for (int i = 0; i < X_B + 2; i++) {
                Vector3 v1 = L_Bound + Vector3(hX * (i - 1), hY * (j - 1), hZ * (k - 1));
                Vector3 v2 = L_Bound + Vector3(hX * i, hY * j, hZ * k);

                cells[k][j][i].curIn = 0;
                cells[k][j][i].stIn = 0;
                cells[k][j][i].o_f = bound[0]->belon(v1, v2) || bound[1]->belon(v1, v2) || bound[2]->belon(v1, v2);
            }
        }
    }
    ///////////////////////////////////////////////////////////////////////////
    for (int q = 0; q < Num_Tube; q++) {
        if (sys[q] != 0) {
            int molR = sys[q]->getRank();
            for (int z = 0; z < molR; z++) {
                Component* b = sys[q]->getComponent(z);
                Vector3 v = b->getPos() - L_Bound;

                int I = v.x / hX + 1;
                int J = v.y / hY + 1;
                int K = v.z / hZ + 1;

                if (I < (X_B + 2) && J < (Y_B + 2) && K < (Z_B + 2) && I > 0 && J > 0 && K > 0) {
                    if (cells[K][J][I].curIn < (int) MAX_PARTICLES) {
                        cells[K][J][I].inCell[cells[K][J][I].curIn] = b;
                        cells[K][J][I].stIn += 1;
                        cells[K][J][I].curIn += 1;
                    } else {
                        fprintf(Log, "Too much particles in cell[%d][%d][%d]: %d (limit is %d)\n", K, J, I, cells[K][J][I].curIn, MAX_PARTICLES);
                        return 0;
                    }
                }
            }
        }
    }
    ///////////////////////////////////////////////////////////////////////////
    float tau = myin.readFloat(TAU, &ok);
    if (!ok) {
        fprintf(Log, "cannot read %s\n", TAU);
        return 2;
    }
    Particle::t = tau;
    Bound::t = tau;
    Flow1::t_step = tau;

    float mRlxTime = ceil(myin.readFloat(RLXT, &ok) / tau);
    if (!ok) {
        fprintf(Log, "cannot read %s\n", RLXT);
        return 2;
    }

    //////////////////////////////////////////////////////

    int OutStep = myin.readFloat(OSTP, &ok);
    if (!ok) {
        fprintf(Log, "cannot read %s\n", OSTP);
        return 2;
    }

    int NUM_OF_ITER = myin.readFloat(NUMITER, &ok);
    if (!ok) {
        fprintf(Log, "cannot read %s\n", NUMITER);
        return 2;
    }

    fclose(InfoFile);

    /////////////////////////////////////////
    FILE* outDynFile = fopen("./IO/Dynamics.dyn", "w");
    if (!outDynFile) {
        fprintf(Log, "Cannot open file ./IO/Dynamics.dyn\n");
        return -2;
    } else {
        fprintf(Log, "File opened: Dynamics.dyn\n");
    }

    fprintf(outDynFile, "<Iterations ScalMin=\"%f\" ScalMax=\"%f\">\r\n", tmk / 2, tmk * 2.5);

    fprintf(Log, "������ �������������.\n");

    Measurer measurer = Measurer(mIn.Info, sub, mIn.Tube, mDscr, L_Bound, R_Bound);
    measurer.outHead();
    measurer.outSub();
    measurer.outChl();

    FILE* newValue = fopen("./IO/Value.vl", "w");
    fprintf(newValue, "Titel:\t%s\r\n", mIn.Info.ModelName);
    fprintf(newValue, "User:\t%s\t%s\r\n\r\n", mIn.Info.Users[0].FName, mIn.Info.Users[0].LName);

    fprintf(newValue, "Substance:\t");
    for (int y = 0; y < sub.NUM_OF_COMPONENT - sub.TubeComp; y++) {
        fprintf(newValue, "%s: %i %%", sub.idName[y].CompName, (int) mDscr[y].Portion);
    }
    fprintf(newValue, "\n");

    fprintf(newValue, "Tube:\t%s \r\n", mIn.Tube->mName);
    fprintf(newValue, "Pressure(Pa): %.2e\r\n", prs);
    fprintf(newValue, "Time(ns):\t%.2e\tTime Step(ns): %.2e\r\n", NUM_OF_ITER * tau, tau);
    fprintf(newValue, "\r\n");
    fprintf(newValue, "-----------------------------------------------------------------------------\r\n");
    fprintf(newValue, "Time(ns)\t\tStream(kgm/s)\tInTemp(K)\tTubeTemp(K)\tOutTemp(K)\t\r\n");
    fflush(newValue);

    Measur measur;
    _ZeroStruct(&measur);

    int oCh = 0;
    int mITER = 0;

    for (int ITER = 0; ITER < NUM_OF_ITER; ITER++) {
        // main cycle, 0 node

        if (ITER % IN_PRD == 0) {
            int Div = inPool.addPool(moleculesCount);
            for (int i = 0; i < Div; ++i) {
                sys[i] = inPool.getAddedMolecula(i);
            }
            moleculesCount += Div;
        }
        ////////////////////////////////////////////////////

        Boxing(Num_Tube, moleculesCount, X_B, Y_B, Z_B, hX, hY, hZ, L_Bound);

        if (ITER % OutStep == 0) {

            printf("ITER Num %i \n", ITER);

            fprintf(outDynFile, "<Iteration Num_of_iter=\"%i\" Time=\"%f\">\r\n", ITER, ITER * tau);
            float Eff_T = 0;
            int outCh = 0;
            for (int k = 0; k < moleculesCount; k++) {
                if (sys[k] != 0) {
                    int R = sys[k]->getRank();
                    for (int l = 0; l < R; l++) {
                        Component* buf = sys[k]->getComponent(l);
                        const Vector3& pos = buf->getPos();
                        const Vector3& vel = buf->getVel();
                        int c_id = buf->getCompID();

                        Eff_T = vel.length_2() * buf->getM() / 3 / 1.38e-14f;
                        fprintf(outDynFile, "ID|%i|CompID|%i|StateID|%i|POS|%f|%f|%f|SCAL|%f|;\r\n", outCh, c_id, 0, pos.x, pos.y, pos.z, Eff_T);
                        outCh++;
                    }
                }
            }
            fprintf(outDynFile, "</Iteration>\r\n");
        }

        // allocations speedup
        int max_rnk = 0;
        for (int k = 0; k < Z_B + 1; k++) {
            for (int j = 0; j < Y_B + 1; j++) {
                for (int i = 0; i < X_B + 1; i++) {
                    int rnk = cells[k][j][i].curIn + cells[k][j][i + 1].curIn + cells[k][j + 1][i + 1].curIn + cells[k + 1][j][i + 1].curIn;
                    if (rnk > max_rnk) {
                        max_rnk = rnk;
                    }
                }
            }
        }

        mpi_toSend * packet = new mpi_toSend[max_rnk];

        for (int k = 0; k < Z_B + 1; k++) {
            for (int j = 0; j < Y_B + 1; j++) {
                for (int i = 0; i < X_B + 1; i++) {
                    int rnk = cells[k][j][i].curIn + cells[k][j][i + 1].curIn + cells[k][j + 1][i + 1].curIn + cells[k + 1][j][i + 1].curIn;
                    if (rnk <= 1) {
                        continue;
                    }

                    Cell* cellBatch[4] = {
                            &cells[k][j][i],
                            &cells[k][j][i + 1],
                            &cells[k][j + 1][i + 1],
                            &cells[k + 1][j][i + 1],
                    };
                    int shift[4] = { 0, 0, 0, 0 };

                    for (unsigned int cellIdx = 0, q = 0; cellIdx < _NElem(cellBatch); ++cellIdx) {
                        const Cell* cell = cellBatch[cellIdx];
                        for (int c = 0; c < cell->curIn; ++c, ++q) {
                            cell->inCell[c]->getPack(packet[q]);
                        }
                        shift[cellIdx] = cell->curIn;
                    }

                    mpiCalc(rnk, shift, packet, tau);

                    for (unsigned int cellIdx = 0, q = 0; cellIdx < _NElem(cellBatch); ++cellIdx) {
                        Cell* cell = cellBatch[cellIdx];
                        for (int c = 0; c < cell->curIn; ++c, ++q) {
                            cell->inCell[c]->addAccel(packet[q].Accel);
                        }
                    }
                }
            }
        }
        delete[] packet;
        ////////////////////////////////////
        inPool.normVel(5);

        for (int k = 0; k < Z_B + 2; k++) {
            for (int j = 0; j < Y_B + 2; j++) {
                for (int i = 0; i < X_B + 2; i++) {
                    if (cells[k][j][i].o_f) {
                        for (int z = cells[k][j][i].stIn; z < cells[k][j][i].curIn; z++) {
                            if (cells[k][j][i].inCell[z] != 0) {
                                bound[0]->bCtrl(cells[k][j][i].inCell[z]);
                                bound[1]->bCtrl(cells[k][j][i].inCell[z]);
                                bound[2]->bCtrl(cells[k][j][i].inCell[z]);
                                bound[3]->bCtrl(cells[k][j][i].inCell[z]);
                            }
                        }
                    }
                }
            }
        }

        float ttemp = 0;
        float otemp = 0;
        float itemp = 0;

        int tch = 0;
        int och = 0;
        int ich = 0;
        float koeff = 2.0 / 3.0 / 1.38e-14;

        for (int k = Num_Tube; k < moleculesCount; k++) {
            if (sys[k] != 0) {
                if (sys[k]->getPos().x > Sht && sys[k]->getPos().x < Sht + mIn.Tube->Len) {
                    ttemp += sys[k]->Ek;
                    tch++;
                }

                if (sys[k]->getPos().x > Sht + mIn.Tube->Len) {
                    otemp += sys[k]->Ek;
                    och++;
                }

                if (sys[k]->getPos().x < Sht - 0.2) {
                    itemp += sys[k]->Ek;
                    ich++;
                }
            }
        }

        if (tch > 10) {
            ttemp = ttemp / tch * koeff;
        }
        if (ich > 10) {
            itemp = itemp / ich * koeff;
        }
        if (och > 10) {
            otemp = otemp / och * koeff;
            oCh++;
        }

        measur.tempOutMA[ITER % MAP] = otemp;
        measur.tempInMA[ITER % MAP] = itemp;
        measur.tempTubeMA[ITER % MAP] = ttemp;
        measur.streamC = inPool.getStream();

        bool gM = tau * ITER > 1000e-5;

        if (gM) {
            measur.tempInAv += measur.tempInMA[ITER % MAP];
            measur.tempTubeAv += measur.tempTubeMA[ITER % MAP];
            measur.tempOutAv += measur.tempOutMA[ITER % MAP];
            measur.streamAv += measur.streamC;

            measur.prsInAv += measur.prsInMA[ITER % MAP];

            mITER++;
        }

        if (ITER > NUM_OF_ITER - 2) {
            fprintf(newValue, "----------------------------------------------------\n");
            fprintf(newValue, "Average:\t\t%.2e\t%.2f\t\t%.2f\t\t%.2f\r\n", measur.streamAv / mITER, measur.tempInAv / mITER, measur.tempTubeAv / mITER, measur.tempOutAv / (oCh + 1));
            fflush(newValue);

            fprintf(outDynFile, "</Iterations>\r\n");
            fclose(Log);
            fclose(outDynFile);
        }

        //////////////////////////////////////
        if (ITER % OutStep == 0 && gM) {
            float outTTemp = getMA(measur.tempTubeMA);
            float outOTemp = getMA(measur.tempOutMA);
            float outITemp = getMA(measur.tempInMA);
            fprintf(newValue, "%.2e\t\t%.2e\t%.2f\t\t%.2f\t\t%.2f\r\n", ITER * tau, measur.streamC, outITemp, outTTemp, outOTemp);
            fflush(newValue);
        }
    }
    printf("All done!\n");
    return 0;
}
